Systems and Methods for Automatic Profiling of Network Event Sequences

ABSTRACT

Systems and methods are disclosed that profile event sequences by creating a mixture model from the event sequences; estimating parameters for the mixture model; and applying the mixture model to profile the event sequences.

BACKGROUND

This invention relates to automatic profiling of network event sequences.

The behavior of network objects, such as flows, sessions, hosts, and end users, can often be described by sequences of communication events in the time domain. Understanding the behavior of networking objects such as traffic flows, sessions, hosts, or users are essential in many applications of network measurement and monitoring. Such behavior can often be described by event sequences, where by event sequence refers to a series of events that i) they are affiliated with the same entity, ii) each event is identified by a symbol, iii) symbols take limited discrete values. While many event sequences can be found in a variety of networking scenarios, three example cases include:

TCP SYN/FIN/RST sequences: The TCP protocol signals the start and the end of a TCP connection with packets that are distinguished by flags in the header. The first packet has a SYN flag set; the last usually has the FIN flag set. A TCP connection can also be terminated by a packet with the RST flag set in the header. The arrival of SYN, FIN and RST in a TCP connection forms an event sequence.

SIP-based VoIP call sessions: SIP (Session Initiation Protocol) is the defactor signaling protocol for VoIP services. Because VoIP relies on SIP to setup and tear down call sessions, each session contains SIP control messages like INVITE, ACK and BYE. Such control messages in each session form an event sequence.

Wi-Fi user sessions: In Wi-Fi wireless networks, a user needs to establish a wireless connection with the nearby access point (AP) to access the Internet. In this case, a user session can be defined as the duration between the user joining and leaving the wireless network. A user session is established, secured and terminated through control message exchange between the user and AP. All such messages within a user session become an event sequence.

As a common problem in network measurement and monitoring tasks, operators often seek simple and effective solutions to understand the diverse behavior hidden in a large amount of event sequences—including both the massive normal behavior and potentially the behavior of a small proportion of anomalies.

The system provides for a complete profiling of massive event sequences. Accurate yet information-compact profiling of such event sequences is critical for many network measurement and monitoring tasks. The system can handle the multi-dimensional behavior exhibited by event sequences. That is, the sequence behavior can not be fully described by a single variable or distribution. Instead, it possesses at least two types of important properties: sequential patterns constituted by symbols, and duration between events. In practice, even one type of properties might be difficult to be described precisely. E.g., although the aforementioned TCP SYN/FIN/RST sequence has only three discrete symbols, the system can handle a large number of patterns with the longest pattern having in excess of 100 symbols (based on a trace collected at a gateway router). For duration related applications, the system can provide a precise profiling. Taking the VoIP sequence in FIG. 1 as an example: the duration between ACK and BYE is the actual call duration. This duration is heavy-tailed and ranges from 0 second to more than 2 hours. When both sequential patterns and duration are concerned, the complexity is conceivably much higher.

There are a large number of prior work on modeling or mining sequence-alike data in areas as diverse as speech recognition, bio-informatics, database, database mining and system. In speech recognition field, applied hidden Markov models (HMM) to cluster a string of acoustic units. In database and database mining field, mining sequential patterns have been studied intensively in the past years. Many of these work focus on discovering rules instead of modeling. For example, Internet users' navigation records are modeled by using a mixture model consisting of first-order Markov chains. In system areas, a relevant research direction is to detect anomalies (or intrusion) via mining sequential system states such as system calls.

In networking areas there are few studies on sequence-alike data analysis and applications. Some recent efforts are devoted to inferring properties of traffic flow. Recently, certain work applies various data mining techniques to identify significant patterns or insignificant anomalies from various traffic entities such as flows, Internet backbone traffic and host communication patterns.

SUMMARY

Systems and methods are disclosed that profile event sequences by creating a mixture model from the event sequences; estimating parameters for the mixture model; and applying the mixture model to profile the event sequences.

Advantages of the preferred embodiment may include one or more of the following. The system supports profiling the behavior hidden in massive event sequences. On a given sequence set, the system automatically learns a mixture model which fully captures the sequence behavior including both event order and duration between events. The learned mixture model is information-compact as it classifies sequences into a set of behavior templates, which are described by multiple continuous-time Markov chains. The model parameters are learned in an iterative procedure based on the Expectation Maximization algorithm. Various network management applications can be done: a visualization tool for network administrators to conduct exploratory traffic analysis, and an efficient anomaly detection mechanism. The accuracy of the system and the efficiency of the two applications have been validated with real data traces collected in three network scenarios: TCP packet traces, VoIP calls, and syslog traces in wireless networks. The system produces a convenient manifestation of the massive behavior within a set of event sequences. The method inherently classifies sequences into so-called elephants and mice categories. This opens a door for a variety of exploratory data analysis tasks. The visualization tool can show the produced compact category information and is valuable for network operators to have an intuitive and easy understanding on behavior patterns in large amount of traffic data. Secondly, because the method provides a way to profile activities of network objects, it can be applied to distinguish normally-behaved objects and abnormally-behaved ones. Such a capability is essential for the design of anomaly detection schemes. The continuous-time Markov chain based mixture model captures both sequential pattern and time information. The method is particularly suitable for networking problems. The system can target network behavior that can be described by a sequence of symbols which take values from a limited number of discrete values. The system provides system operators with tools to conveniently monitor massive traffic.

BRIEF DESCRIPTION OF THE DRAWINGS

FIG. 1 shows three exemplary event sequence cases.

FIG. 2 shows an exemplary process to model event sequences.

FIG. 3 shows an exemplary continuous time Markov chain for TCP SYN/FIN/RST sequences.

FIG. 4 shows an exemplary EM algorithm.

FIG. 5 shows a mixture model of two CTMCs and an observed event sequences.

FIG. 6 compares the two histograms on TCP SYN/FIN/RST sequences.

FIG. 7-8 show exemplary histograms for VoIP call sessions.

FIG. 9 shows a comparison of histogram for discrete symbols.

FIG. 10 shows an exemplary visualization of CTMCs.

FIG. 11 shows an anomaly score chart for TCP sequences.

FIG. 12 shows a visualization of significant CTMCs for a period of time.

FIG. 13 shows a second embodiment of a learning model.

DESCRIPTION

FIG. 2 shows an exemplary process for profiling network event sequences. The procedure of event sequence modeling first creates a mixture model (10). In one embodiment, the model includes multiple first-order Markov chains and only transition probabilities are modeled by parameters. Next, the process estimates one or more parameters using the Expectation Maximization process (20). The mixture model is then deployed in an application (30). In one application, visualization of sequences is done (40). The sequences assigned to each Markov chain is visualized by considering the transition probabilities which can be color coded. In another application, the system uses the profile to model the behavior of network objects (50).

The method of FIG. 2 can capture the behavior of massive sequences with compact information. Given a large sequence dataset, the method automatically learns a mixture model which consists of multiple continuous-time Markov chains. The model is general in the sense that it requires little knowledge of the targeting sequences except for the number of discrete symbols. The model is expressive because its structure enables it to capture simultaneously the two aspects of sequence behavior, i.e., sequential symbol pattern and duration between symbols. Also, the model has flexibility of adjusting its sensitivity on either aspect to match application requirements. In model learning process, the system uses maximum likelihood estimation to find values for the model parameters. The estimation procedure is developed from the Expectation-Maximization (EM) algorithm. Although the method does not require prior domain knowledge of the targeting sequences, such knowledge can be utilized to improve model accuracy and to reduce computational complexity. The method is evaluated by applying it to the three aforementioned sequence examples. Evaluation experiments are based on three real data traces, including a TCP packet trace, a syslog trace in wireless networks, and VoIP calls logged at a VoIP testbed set up on PlanetLab. Experiment results demonstrate that the method produces accurate profiles for the real data.

Next, the mixture model approach is discussed in more detail. In a continuous-time Markov chain for TCP SYN/FIN/RST sequences, each transition is labeled by its probability and rate. The rate for any transition adjacent to the “start” and “end” states are infinite because such transitions are considered to take zero time.

The system learns a model from input sequence data. Once the model is learned, it is expected to mimic the sequence behavior. The basic model assumes that an event sequence is independently generated from one out of C Continuous Time Markov Chains (CTMCs), and the sequence is mapped onto a path in it. The approach is also called mixture-model based since the model contains multiple components. Each component characterizes one type of sequence behavior so that the entire model is able to characterize the diverse behavior in massive sequences.

Turning now to the structure of CTMCs, the CTMCs in the mixture model have the same structure. In each CTMC, there is a corresponding state for every discrete symbol value appearing in sequences. In addition, CTMC includes two special states: “start” and “end”. CTMC is full meshed, i.e., any state can transit to another state, including the state itself. The only exceptions are: the “start” state has no incoming transitions and the “end” state has no outgoing transitions. For each transition, both its probability and duration are modeled. For the transition from state u and v, its probability is denoted by q_(uv); its delay is assumed to follow an Exponential distribution with a parameter λ_(uv), which is also called transition rate. As an example, FIG. 3 shows a CTMC for the aforementioned TCP SYN/FIN/RST sequence.

Given the above structure, any event sequence can be mapped to a path from the “start” state to the “end” state on CTMC. Also, creating such a CTMC only requires to know the total number of discrete symbols. These two features make the model applicable for various event sequences.

To formally describe the problem, S denotes the input sequence set. s denote one sequence instance in S. N denotes the size of S. A complete description of s is a tuple {b₁(s), t_(1,2)(s), b₂(s), t_(2,3)(s), b₃(s), . . . , b_(l(s))(s)}, where b_(i)(s) is the i-th symbol in s, t_(i,i+1)(s) is the duration between the i-th and i+1-th symbol, and l(s) is the sequence length. M denotes the number of states in each CTMC. M includes the “start” and the “end” state. Θ denotes the entire parameter set of the mixture model. θ^(c) denotes parameters in CTMC c (c=1, . . . , C). θ^(c) includes all those transition parameters q_(uv) and λ_(uv). If necessary, a superscript is added to such transition parameters to specify which CTMC they belong to. For example, q_(uv) ^(c) denotes the transition probability from state u to v in CTMC c.

Given the above notations, the probability of observing a sequence s is given by

${p\left( s \middle| \Theta \right)} = {\sum\limits_{c = 1}^{C}{\pi_{c}{p\left( s \middle| \theta^{c} \right)}}}$

where π_(c) is called mixture probability. It is also a parameter included by Θ. π_(c) shall be interpreted as the a priori probability of assigning a sequence to CTMC c. π_(c) must satisfy the constraint

${\sum\limits_{c = 1}^{C}\pi_{c}} = 1.$

Initially, all parameters in the model, Θ, are unknown. The goal is to train the model with input data and find those parameters that best explain sequence behavior in the data. In the next section, an exemplary algorithm for parameter estimation is discussed.

The Expectation Maximization (EM) algorithm is used to derive parameters in the mixture model. EM algorithm is a general method for finding the maximum likelihood estimate for model parameters when the observed data is incomplete. The procedure of EM algorithm is shown in FIG. 4 and presented below

1. Let l = 0 and initialize Θ_(l) 2. do 3. ${Q\left( {\Theta,\Theta_{l}} \right)} = {\sum\limits_{m \in \partial}\; {\log \mspace{11mu} p\; \left( {S,{m\Theta}} \right)\mspace{11mu} p\mspace{11mu} \left( {m\Theta_{l}} \right)}}$ 4. Compute Θ_(l+1) = argmax_(Θ)Q(Θ, Θ_(l)) 5. l = l + 1 6. until (convergence criterion is satisfied)

In the above algorithm, the estimate for Θ is updated in an iteration manner. In the l+1-th loop, a new estimate for Θ is obtained and denoted by Θ_(l+1). Θ_(l+1) is computed based on training data and Θ_(l), which is the estimate obtained in the l-th loop. The iteration procedure is guaranteed to converge to a solution Θ that yields a local maximum value for the target likelihood function log(Θ|S). In applying EM algorithm, the system determines Q(Θ,Θ_(l)) and maximizes Q(Θ,Θ_(l)) in each loop. These are the so-called expectation step and maximization step. The following equations hold for a mixture model

$\begin{matrix} {{Q\left( {\Theta,\Theta_{l}} \right)} = {{\sum\limits_{c = 1}^{C}{\sum\limits_{s \in S}{{\log \left\lbrack {p_{c}\left( s \middle| \Theta \right)} \right\rbrack}{p\left( {\left. c \middle| s \right.,\Theta_{l}} \right)}}}} + {\sum\limits_{c = 1}^{C}{\sum\limits_{s \in S}{{\log \left( \pi_{c} \right)}{p\left( {\left. c \middle| s \right.,\Theta_{l}} \right)}}}}}} & {(1)(2)} \end{matrix}$

where p(c|s,Θ_(l)) is the posterior probability of assigning s to CTMC c. It can be computed via Bayes' theorem that

${p\left( {\left. c \middle| s \right.,\Theta_{l}} \right)} = \frac{\pi_{c}{p\left( s \middle| \theta^{c} \right)}}{\sum\limits_{k = 1}^{C}{\pi_{k}{p\left( s \middle| \theta^{k} \right)}}}$

wherein and Θ^(k) are π_(k) known from Θ_(l), so p(c|s,Θ_(l)) is determinable.

In ((2)), p_(c)(s|Θ) is the probability of generating a sequence s from CTMC c which takes parameters from Θ. Its mathematical form can be derived is as following

$\begin{matrix} {{{p_{c}\left( s \middle| \Theta \right)} = {{p_{c}\left( s \middle| \theta^{c} \right)} = {\prod\limits_{j = 1}^{{l{(s)}} - 1}q_{s_{j}}^{c}}}},{s_{j + 1}{p\left( {{\lambda_{b_{j}}^{c}(s)},{b_{j + 1}(s)},{t_{j,{j + 1}}(s)}} \right)}}} & (3) \end{matrix}$

where

p(λ_(b_(j))^(C)(s), b_(j + 1)(s), t_(j, j + 1)(s))

is the a priori probability of observing a transition delay t_(j,j+1)(s) given that the delay follows an Exponential distribution with a parameter

λ_(b_(j))^(C)(s), b_(j + 1)(s).

Such a probability is determined by the following function

$\begin{matrix} {{p\left( {\lambda,t} \right)} = \left\{ \begin{matrix} {{^{{- \lambda}\; \frac{t}{\Delta \; t}}\left( {1 - ^{- \lambda}} \right)},} & {t \leq t_{\max}} \\ {{^{{- \lambda}\; \frac{t_{\max}}{\Delta \; t}}\left( {1 - ^{- \lambda}} \right)},} & {t > t_{\max}} \end{matrix} \right.} & (4) \end{matrix}$

where Δt and t_(max) are two user-defined constants. p(λ,t) returns a fixed yet small probability when the transition delay t is larger than t_(max), When t is less than or equal to t_(max), p(λ,t) returns the accumulated Exponential density probability over the time interval

$\left\lbrack {\frac{t}{\Delta \; t},{\frac{t}{\Delta \; t} + 1}} \right\rbrack.$

In fact, alternative function forms can be chosen for p(λ,t). Here ((4)) is chosen is because it is more tractable. The parameters t_(max) and Δt should be carefully chosen because they affect the model accuracy and its sensitivity with respect to different properties of sequences.

Now Q(Θ,Θ_(l)) in ((1) is maximized. The two parts on the right side are independent and can be maximized separately. To maximize the first part, by substituting ((3)) ((4)) into the first part:

$\begin{matrix} {\sum\limits_{c = 1}^{C}\; {\sum\limits_{s \in S}\; {\sum\limits_{j = 1}^{{l{(s)}} - 1}\; {\left( {{\log \; {q_{b_{j}}^{c}(s)}},{{b_{j + 1}(s)} - {\lambda_{b_{j}}^{c}(s)}},{{{b_{j + 1}(s)}{t_{j,{j + 1}}^{*}(s)}} + {\log \; {\lambda_{b_{j}}^{c}(s)}}},{b_{j + 1}(s)}} \right){p\left( {\left. c \middle| s \right.,\Theta_{l}} \right)}}}}} & (5) \end{matrix}$

where

${t_{j,{j + 1}}^{*}(s)} = \frac{t_{j,{j + 1}}(s)}{\Delta \; t}$

if t_(j,j+1)(s)≦t_(max); otherwise

${t_{j,{j + 1}}^{*}(s)} = {\frac{t_{\max}}{\Delta \; t}.}$

Note that the derivation of ((5)) uses an approximation 1−e^(−λ)≈λ.

To maximize ((5)) with respect to each transition parameter, a transition from state u to state v is done in CTMC c. The transition may appear in a sequence s for an arbitrary number of times. I(s, u, v) denotes the frequency that the transition appears in s. Γ(s, u, v) denotes the aggregate duration for all such transitions. Again, in computing F(s, u, v) the previously defined t_(u,v)*(s) is used whenever t_(u,v)(s) is involved. Since transition probabilities are constrained by

${\sum\limits_{v = 1}^{M}\; q_{uv}^{c}} = 1$

for every u, the Lagrange multiplier method is used to maximize ((5)) with respect to q_(u,v) ^(c) and λ_(u,v) ^(c). The following equations are obtained:

$\begin{matrix} {{q_{u,v}^{c} = \frac{\sum\limits_{s \in S}\; {{I\left( {s,u,v} \right)}{p\left( {\left. c \middle| s \right.,\Theta_{l}} \right)}}}{\sum\limits_{i = 1}^{M}\; {\sum\limits_{s \in S}\; {{I\left( {s,u,i} \right)}{p\left( {\left. c \middle| s \right.,\Theta_{l}} \right)}}}}}{\lambda_{u,v}^{c} = \frac{\sum\limits_{s \in S}\; {{I\left( {s,u,v} \right)}{p\left( {\left. c \middle| s \right.,\Theta_{l}} \right)}}}{\sum\limits_{i = 1}^{M}\; {\sum\limits_{s \in S}\; {{\Gamma \left( {s,u,i} \right)}{p\left( {\left. c \middle| s \right.,\Theta_{l}} \right)}}}}}} & {(6)(7)} \end{matrix}$

Next, the second part of the right side in ((1)) is maximized with respect to π_(c). Again the Lagrange multiplier method is applied because of the constraint

${\sum\limits_{c = 1}^{C}\; \pi_{c}} = 1.$

π_(c) is solved as

$\begin{matrix} {\pi_{c} = \frac{\sum\limits_{s \in S}\; {p\left( {\left. c \middle| s \right.,\Theta_{l}} \right)}}{N}} & (8) \end{matrix}$

The three equations, ((6)) ((7)) and ((8)), comprise the parameter estimation process at the l+1-th loop of the iteration procedure.

Using protocol knowledge to fix certain parameters, certain parameters are fixed to constants by using the prior knowledge on involved networking protocols. This not only improves the model accuracy, but also reduces computational complexity because the fixed parameters no longer participate in the parameter estimation procedure. Both transition probabilities and rates can be fixed. Below are two commonly used strategies:

if a transition from state u to state v is impossible in reality, the corresponding transition probability can be set to zero; For example, in the TCP case, it is known that a FIN packet is never followed by a SYN in the same connection. Accordingly, the probability for the transition from “FIN” state to “SYN” state is set to zero. Clearly, in this case the rate for the involved transition also needs to be fixed.

if a transition from state u to state v is always finished in very short time, the corresponding transition rate can be fixed to a large constant; for the VoIP example: RING and OK messages are always sent together by the VoIP server. Thus their corresponding transition rate is set to a large constant.

Because the model is a collection of multiple CTMCs, prior to any probability (or rate) fixing for the transition between state u and v, the number of CTMCs this operation should be applied to needs to be selected. In principle, such a decision depends on the confidence level on the domain knowledge. A high confidence level implies that the knowledge always holds in reality and no exceptions will happen; if this is true, the parameter fixing operation should be applied to all the CTMCs. Otherwise, only a subset of CTMCs can have their parameters fixed in the hope that the other CTMCs can still adjust their parameters to model exceptions.

Next, the selection of initial values for free parameters will be discussed. After certain parameters are fixed with domain knowledge, the rest are undecided parameters and will be estimated via the iterative procedure. Prior to the iteration procedure, initial values are assigned to them. In practice, random numbers are used as the initial values. A normalization operation is added to ensure that those initial values satisfy constraints such as

${\sum\limits_{c = 1}^{C}\; \pi_{c}} = 1$ and ${\sum\limits_{j = 1}^{M}\; q_{ij}^{c}} = 1.$

As EM algorithm only guarantees local optimality, accuracy of the method is affected by the initial values. Accordingly, the model training process is repeated multiple times (10 in this work) with different random seeds. The model having the best compliance with the convergence criterion (see Line 6 in FIG. 4). Regarding the convergence criterion, Q(Θ,Θ_(l)) is computed upon the finish of the l-th loop. If

$\frac{{{Q\left( {\Theta,\Theta_{l}} \right)} - {Q\left( {\Theta,\Theta_{l - 1}} \right)}}}{Q\left( {\Theta,\Theta_{l - 1}} \right)} < \delta$

(δ is a user-defined threshold which is set to 0.05% in the experiments), the iteration is terminated.

Next, the selection of C is done. C (the total number of CTMCs) is selected prior to the iteration procedure. Selecting an appropriate value for C is needed because it balances model accuracy with complexity. The method uses the Akaike Information Criterion (AIC) which is defined as

AIC(Θ)=−2 log p(S|Θ)+2U

where U is the free parameter number in the model. It is easy to show that the total parameter count of the model is C(2M²−7M+6)−1. After using protocol knowledge to fix certain parameters, U should be much smaller than the total parameter count. The introduced AIC measures the goodness of fit for models. A chosen C that minimizes AIC(θ) is interpreted as a good balance between model accuracy and complexity. In applying AIC, an upper bound for C is selected and then the system searches below the upper bound for a choice of C that minimizes AIC(O). How to choose the upper bound depends on the number of states in each CTMC, as well as the prior knowledge on how heterogeneous the sequence behavior could be. A high upper bound should be chosen if CTMC has many states or the sequence behavior are highly heterogeneous.

The major computational complexity is caused by the three equations ((6)) ((7)) and ((8)). Their required complexity is to multiply N by the total number of free parameters, which is bounded by O(CM²). Prior to solving these three equations, ((3)) and ((4)) must be computed. While Equation ((3)) requires O(CN), Equation ((4)) involves an Exponential function which can be replaced by a lookup table in implementation. So the complexity for ((4)) is O(M²). Overall, the computational complexity in one iteration is upper bounded by O(NCM²+CN+M²)≈O(NCM²). This complexity is multiplied by the number of iterations. This number is less than twenty in all the experiments.

As mentioned above, the presented model is expected to characterize two aspects of the sequence behavior: the sequential pattern constituted by symbols, and durations between symbols. It is important to understand the model sensitivity with respect to each aspect because this decides the capability of the model in distinguishing sequences with different behavior. In this section the model sensitivity is analyzed. An example shows how the model sensitivity is affected by parameter selection. The model sensitivity can be tailored according to application requirements.

FIG. 5 shows a mixture model of two CTMCs and an observed event sequence s. The two CTMCs are assumed to have equal mixture probability. In the two CTMCs only those edges with non-zero transition probabilities are depicted. Each edge is labeled by its transition probability and rate. The likelihood for s being generated by CTMC 1 and 2, respectively, is determined. In CTMC 1, the expected transition delay from state 1 to itself is 1, and the expected delay from state 1 to 2 is 1/0.1=10. These two numbers are consistent with their respective observations in s. Therefore, if transition delay is the primary concern, s is more likely generated by CTMC 1 since in CTMC 1 the property of transition delay fits observations well. On the other hand, if the order of symbols is the primary concern, s has a higher chance to be generated from CTMC 2. Because in CTMC 2 the self transition probability at state 1 is much higher (0.9) than in CTMC 1, which better explains the observed multiple consecutive symbol “1”. As illustrated by this example, the model has a trade-off between being more sensitive to the order of symbols or to transition durations.

The model sensitivity can be controlled by two parameters: t_(max) and Δt (see ((4))). From ((4)), it is clear that the difference of transition delay only matters when the duration is less than or equal to t_(max). If t_(max) is small, the model tends to be insensitive to variance of transition delay. In the last example, if t_(max)=100, which is larger than any transition delay in s, and Δt=1, then the probability of s being generated from each CTMC is as follows

p(s|θ ¹)=[0.1×(e ⁰ −e ⁻¹)]²×0.9×(e ⁻¹ −e ^(−1.1))

p(s|θ ²)=[0.9×(e ⁰ −e ⁻¹)]²×0.1×(e ⁻¹⁰ −e ⁻¹¹)

There is p(s|θ¹)?p(s|θ²). Since the two CTMCs have equal mixture probabilities, the model has a higher posterior probability to assign s to CTMC 1 than to CTMC 2. On the other hand, if Δt=10 and t_(max) is unchanged, there is p(s|θ¹)=p(s|θ²), the assignment becomes different. This example gives the following insight: while a smaller Δt makes the model more sensitive to transition delay, a larger Δt makes it more sensitive to the order of symbols. In an extreme scenario, to ascertain whether the order of symbols follows certain protocol logic, then t_(max)=0. In the other extreme scenario where the system is only interested in duration, t_(max) is made larger than any possible duration values and Δt is close to zero.

Lastly, it is possible to specify different t_(max) and Δt for different transitions. This feature is particularly useful when transitions have unequal importance. E.g., if a transition is not within the concern, its t_(max) and Δt are adjusted to let the mixture model ignore the impact of its variance.

The experiments test the three types of event sequences as mentioned above: TCP SYN/FIN/RST sequences, SIP-based VoIP calls and Wi-Fi user sessions. The sample data is extracted from three network traces, respectively. These include a TCP packet trace collected from a gateway router, a VoIP call traces, and a syslog trace in Wi-Fi wireless networks.

The TCP trace is publicly available and is collected from a border router at UCLA during August 2001. The trace is about 16 hours long. From the trace the SYN, FIN, RST packets are extracted and grouped into event sequences by their flow identities, i.e., the tuple (source IP, destination IP, source port, destination port). In total 505K TCP SYN/FIN/RST sequences were extracted. The VoIP trace is collected from a VoIP testbed. An open-source software SIPp for VoIP call generators was used and another open-source software SER was used as the SIP proxy. SIPp was deployed on twenty PlanetLab nodes and generated calls to the SIP proxy. tcpdump was used to capture packets at the proxy and extract all the call sessions. To make the testbed more realistic, intentionally created packet loss and delayed reply were simulated at both call sender and receiver sides. The trace was about 1.5 h the long with 16K call sessions logged. The third trace was syslog data collected from the campus-wide wireless networks in Dartmouth College. The syslog trace contained all the syslog messages exchanged between Wi-Fi wireless cards and their associated access points. In the experiment, only the trace for one access point was used in the first two weeks of 2003. This selected trace includes 2.2M syslog events which are further grouped into 122K user sessions by using the MAC address.

The method was tested against the three traces on a Intel Xeon(R) 1.60 GHz processor. Table 1 summarizes the experiment settings and basic results. The column “Significant CTMCs” is particularly interesting. It reveals how many popular behavior templates the input sequences have. Such significant CTMCs are closely related to the network protocol underlying the data.

TABLE 1 Experiment settings and basic results (“Significant CTMCs”: the number of CTMCs with mixture prob. larger than 5%) Data Δt t_(max) C Significant CTMCs Run-time TCP 1 sec 10 min 6 4  7 min VoIP 1 sec 10 min 12 2 90 sec Wi-Fi syslog 30 sec   2 hour 26 5  8 min

Two metrics were used to evaluate the accuracy of the method: histogram of sequence size and histogram of sequence duration. Both metrics are self-evident. For either metric, A consistency reported by the comparison of the data and fitted model indicates that the model does capture the input sequence behavior. While computing the histograms for the input data is straightforward, it is difficult to analytically derive the histograms from the created models so Monte Carlo simulations was applied. The comparison is conducted for all the three datasets and the results are provided in FIGS. 6-8. All the figures used a logarithmic scale for the sequence duration axis.

FIG. 6 compares the two histograms on TCP SYN/FIN/RST sequences. The sequence size histogram has two spikes at size 1 and 4 (see the left figure). The spike at size 1 is caused by the fact that many TCP connections have a single SYN without any FIN or RST. The spike at size 4 is caused by the most common case wherein a TCP connection has two SYNs and two FINs. For these two spikes, the model seems to give an accurate explanation. The right side of FIG. 6 shows the histogram of sequence duration. Both the real data and fitted model show that more than 70% TCP connections have durations within one second. FIG. 7 plots the histograms for VoIP call sessions. The VoIP call size histogram has a sharp spike at size 6. This corresponds to the most typical VoIP call pattern shown in FIG. 1. On the other hand, the VoIP call duration have a relatively flat distribution. In each case, the model displays similar statistics to the real data. Last, the histogram of Wi-Fi user session size is relatively flat, yet the histogram of session duration is uneven.

The histogram of sequence size is further broken down to examine the histogram of each discrete symbol. The comparison results are shown in FIG. 9. FIG. 9 demonstrates that the model matches the real data extremely well in terms of histograms of discrete symbol values.

Next the applications of the method are discussed with two showcases. The first showcase is a visualization tool that allows network operators to conduct exploratory data analysis; the second showcase is a novel anomaly detection mechanism.

In the exploratory data analysis application, by nature of the mixture model, the method classifies input event sequences into different CTMCs. Each CTMC is a cluster that profiles a group of sequences with similar behavior. This feature is attractive to various exploratory data analysis tasks such as monitoring dominant traffic trend, detecting outliers, developing alternative parsimonious traffic model.

As a showcase, a tool was developed to visualize the clusters produced by the method. Its purpose is to use graphic display to aid network operators' understanding on large volume of event sequence dataset. FIG. 10 shows a visualization of the most significant CTMCs after building models for the TCP SYN/FIN/RST dataset. Transition probabilities in each CTMC are encoded by thickness of the edges. A thicker edge indicates a higher probability and vice versa. Edges with transition probabilities smaller than 5% are pruned in the figure. Transition rates are encoded by the ratio of dash length and dash gap of the edges. A more scattered dashed line indicates a smaller transition rate, i.e., a longer expected transition delay.

FIG. 10 shows the tool's display of the significant clusters after modeling the TCP dataset. The tool has one unique feature: it encodes the transition delay and rate parameters by thickness and scattering of the corresponding edge respectively. This feature, along with the mixture probability numbers appearing on the display, enables the operators to have a coarse, yet immediate estimate of the status of the overall data. For example, by looking at CTMC 1 in the figure, 42% of TCP connections either has only a SYN or has a SYN and a RST. For CTMC 2, 35% of TCP connections have both SYN and FIN. Since the SYN and FIN have non-negligible self transition probabilities in this case, these 35% TCP connections may have multiple consecutive SYNs (FINs). Moreover, if consecutive SYNs occur, their intermediate duration is quite short. While for consecutive FINs, the duration is relatively long. CTMC 4 reveals that 13% of TCP connections have a single or multiple consecutive SYNs. CTMC 4 reveals that 6% of connections have a high chance of having multiple consecutive SYNs, followed by zero or one FIN.

In the above scenario, the tool displays only the significant CTMCs by pruning those insignificant ones. This enables operators to promptly catch the dominant trend hidden in the data and take actions if necessary. For example, it is well known that most TCP RSTs are caused by the implementation of certain Web browsers, network outages, or malicious attacks. If the visualization tool reveals that the popularity of TCP RSTs is increasing sharply, operators are alarmed and they should examine whether a TCP reset attack is launched. On the other hand, the tool can also display those insignificant, yet possibly unusual behavior. Enabling the option or not is purely subject to operators' requirements.

Turning now to anomaly detection, because the modeling approach provides one way to describe the behavior of network objects, it is readily applicable to anomaly detection whose duty is essential to distinguish abnormally-behaved objects from normally-behaved ones. There could be many ways to design an anomaly detection method. Nevertheless, a simple method can use the modeling method and the visualization tool mentioned earlier.

In anomaly detection two sets of event sequences are analyzed: S and Ŝ. While S are known to be clean without anomalies, Ŝ are the test dataset may have anomalies. In reality, Ŝ could be generated by one host or user, or the data collected in a monitoring period. The goal of the anomaly detection technique is to determine whether the behavior of Ŝ is different from S based on certain criteria. If different, S is reported to contain anomalies.

First, a model θ is built for S. θ is considered to represent the normal behavior. In principle another model for Ŝ can be built and compared to Θ. If the two models have significantly different parameter values, Ŝ contains anomalies. Despite its seemingly simplicity, this approach suffers from one drawback: training two models are time-consuming and Ŝ is usually small and insufficient for training a reliable model. Therefore, after obtaining θ, the system extracts the mixture probabilities π_(c) from Θ. π_(c) are the a priori probability for a sequence coming from CTMC c. Such a distribution is denoted by P(CTMC c)=π_(c). This distribution represents the normal behavior. Next, Θ is applied to calculate the membership probability p(c|s,Θ) for every s in Ŝ. p(c|s,Θ) is the probability for assigning s to CTMC c. It satisfies the constraint

${\sum\limits_{c = 1}^{C}\; {p\left( {\left. c \middle| s \right.,\Theta} \right)}} = 1.$

The system computes

${E_{c} = {\sum\limits_{s \in \hat{S}}\; {p\left( {\left. c \middle| s \right.,\Theta} \right)}}},{c = 1},2,\ldots \mspace{11mu},C$

where E_(c) shall be interpreted as the number of sequences in Ŝ that are assigned to CTMC c. It is also called the observed frequency count for CTMC c.

Clearly, if the sequence behavior in Ŝ and in S are similar, the observed frequency counts E₁, . . . , E_(C) should follow the a priori distribution P(CTMC c)=π_(c). Thus this fact is the basis for anomaly detection. Specifically, the system applies χ² (chi-square) statistics to measure the fitness between the observed frequency counts and the distribution, that is, compute

$\chi^{2} = {\sum\limits_{c = 1}^{C}\; \frac{\left( {{{\hat{S}}\pi_{c}} - E_{c}} \right)^{2}}{E_{c}}}$

where |Ŝ| is the size of Ŝ, |Ŝ|π_(c) is the expected frequency count for CTMC c. Because the value of χ² is affected by |Ŝ|, it is normalized in the following way

$C_{c} = \sqrt{\frac{\chi^{2}}{\chi^{2} + {\hat{S}}}}$

The above C_(c) is usually called contingency coefficient. Since it is normalized to be between zero and one, it is readily used to measure how likely an anomaly is detected. It equals zero when Ŝ exhibit identical behavior as S; and it equals one when their behavior do not resemble at all. In practice, a high value of C_(c) indicates Ŝ might be anomalies and require a closer examination.

Overall, the presented technique requires to build a mixture model for a clean set of sequences. This only needs to be done occasionally under a reasonable assumption that the normal behavior does not change frequently. For each testing sample set Ŝ, the technique only needs to compute the membership probability for each sample. The technique is light-weighted and can be performed either online or offline.

For the purpose of evaluation, the technique conducts an offline anomaly detection on the aforementioned TCP dataset. Since there is no clean dataset, the mixture probability distribution computed from the entire dataset is used to represent normal behavior. This is valid as long as the proportion of anomalies in the entire dataset is small. The system then divides the 16-hour trace into 30-minute intervals and computes the contingency coefficient for every interval. The results are shown in FIG. 12 which shows a contingency coefficient larger than 0.4. This value is relatively high and indicates a suspicious interval. The modeling of sequences in the suspicious interval is done so to reveal more information. A model is trained by using the 13K TCP sequences in the suspicious interval. The visualization tool is used to display the modeling results which is shown in FIG. 12. It is interesting to make a visual comparison between FIG. 12 and FIG. 10, which is assumed to reflect the normal behavior. By reading the mixture probability numbers and estimating the thickness of edges on the two figures, in the suspicious interval a much higher proportion of TCP connections contains only a single SYN without FIN or RST. Specifically, in the suspicious interval this percentage is about 60%. While in the normal behavior model the percentage is about 33%.

FIG. 13 shows a second embodiment of a learning model. In this model, the system defines a structure for the mixture model (310) and then estimate model parameters such as by EM (320). In (310), the system determines the number of states in each Markov chain (MC) (312) and then chooses K, the number of MCs in the mixture model (314). In (320), the system initializes parameters with a random guess (322). Next, the process enumerates all possible flows to the MCs and then computes an expect value of log-likelihood (324). The process then updates the random guess by maximizing the expected value of log likelihood (326). If the update value is below a predetermined threshold, the system exits. Otherwise, the process loops to (322) to continue the estimation of the model parameters.

Automatic profiling of massive event sequences is critical for a variety of network measurement and monitoring tasks. The above embodiment automatically captures the behavior hidden in massive event sequences. The approach is based on a mixture model—a collection of continuous-time Markov chains. The most salient feature of the model is that it simultaneously captures both the order of events and duration between events. Prior domain knowledge on the event sequences can be seamlessly integrated into the model to improve accuracy and to reduce complexity. To estimate parameters of the model, an iterative algorithm based on the Expectation Maximization algorithm is used. The method is then applied to multiple network traces, including a TCP packet trace, a VoIP call collection and Wi-Fi syslog and the experimental results demonstrate that the method yields a high consistency with the real data. Furthermore, the method can be applied to two exemplary network monitoring tasks: 1) the visualization tool can be used for conducting exploratory traffic analysis and 2) an anomaly detection scheme that is light-weighted and can be used online or offline.

The invention may be implemented in hardware, firmware or software, or a combination of the three. Preferably the invention is implemented in a computer program executed on a programmable computer having a processor, a data storage system, volatile and non-volatile memory and/or storage elements, at least one input device and at least one output device.

By way of example, a block diagram of a computer to support the system is discussed next. The computer preferably includes a processor, random access memory (RAM), a program memory (preferably a writable read-only memory (ROM) such as a flash ROM) and an input/output (I/O) controller coupled by a CPU bus. The computer may optionally include a hard drive controller which is coupled to a hard disk and CPU bus. Hard disk may be used for storing application programs, such as the present invention, and data. Alternatively, application programs may be stored in RAM or ROM. I/O controller is coupled by means of an I/O bus to an I/O interface. I/O interface receives and transmits data in analog or digital form over communication links such as a serial link, local area network, wireless link, and parallel link. Optionally, a display, a keyboard and a pointing device (mouse) may also be connected to I/O bus. Alternatively, separate connections (separate buses) may be used for I/O interface, display, keyboard and pointing device. Programmable processing system may be preprogrammed or it may be programmed (and reprogrammed) by downloading a program from another source (e.g., a floppy disk, CD-ROM, or another computer).

Each computer program is tangibly stored in a machine-readable storage media or device (e.g., program memory or magnetic disk) readable by a general or special purpose programmable computer, for configuring and controlling operation of a computer when the storage media or device is read by the computer to perform the procedures described herein. The inventive system may also be considered to be embodied in a computer-readable storage medium, configured with a computer program, where the storage medium so configured causes a computer to operate in a specific and predefined manner to perform the functions described herein.

The invention has been described herein in considerable detail in order to comply with the patent Statutes and to provide those skilled in the art with the information needed to apply the novel principles and to construct and use such specialized components as are required. However, it is to be understood that the invention can be carried out by specifically different equipment and devices, and that various modifications, both as to the equipment details and operating procedures, can be accomplished without departing from the scope of the invention itself. 

1. A method to profile event sequences, comprising: a. creating a mixture model that simultaneously characterizes event pattern and time information; b. estimating parameters for the mixture model; and c. applying the mixture model to profile event sequences in computer networks.
 2. The method of claim 1, wherein the mixture model comprises multiple first order Markov chains whose transitions have two parameters: transition probability and transition rate.
 3. The method of claim 1, comprising adding constraints to the mixture model by using domain knowledge.
 4. The method of claim 1, comprising selecting a number of Markov chains in the mixture model by using domain knowledge.
 5. The method of claim 1, comprising estimating parameters using an Expectation Maximization process.
 6. The method of claim 1, comprising visualizing the mixture model.
 7. The method of claim 6, comprising using dotted line to display transition and setting the ratio of dash length and dash gap in proportion to a function of average transition delay.
 8. The method of claim 6, comprising setting the thickness of each transition line in proportion to transition probability
 9. The method of claim 1, comprising detecting network anomalies based on the mixture model.
 10. The method of claim 1, wherein the mixture model characterizes at least two behaviour aspects of massive event sequences.
 11. A system to profile event sequences, comprising: a. a mixture model that simultaneously characterizes event pattern and time information; b. a parameter estimator coupled to the mixture model; and c. a user interface module to apply the mixture model to profile event sequences in computer networks.
 12. The system of claim 11, wherein the mixture model comprises multiple first order Markov chains whose transitions have two parameters: transition probability and transition rate.
 13. The system of claim 11, comprising adding constraints to the mixture model by using domain knowledge.
 14. The system of claim 11, wherein the mixture model comprises a number of Markov chains selected by using domain knowledge.
 15. The system of claim 11, wherein the parameter estimator comprises an Expectation Maximization process.
 16. The system of claim 11, comprising a graphical user interface (GUI) to visualize the mixture model.
 17. The system of claim 16, wherein the GUI uses dotted line to display transition and sets the ratio of dash length and dash gap in proportion to a function of average transition delay.
 18. The system of claim 16, wherein the GUI sets the thickness of each transition line in proportion to transition probability.
 19. The system of claim 11, comprising a verifier to detect network anomalies based on the mixture model.
 20. The system of claim 11, wherein the mixture model characterizes at least two behaviour aspects of massive event sequences. 